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DIRECT-SEARCH  TECH¬ 
NIQUES  DETERMINE  THE 
INDEPENDENT  VARIABLES 
THAT  MINIMIZE  OR 
MAXIMIZE  A  GIVEN 
OBJECTIVE  FUNCTION 


TPBVP  OCCURS  IN  MANY 
PHYSICAL  SCIENCES  - 
HAS  BEEN  EMPHASIZED 
IN  OPTIMAL  CONTROL 
THEORY 


Mathematical  programming  concepts  and  the  study  of  optimization  have  resulted 
in  the  development  of  direct -search  algorithms  that  are  applicable  to  numerous  engineering 
design  problems.  These  techniques  determine  (he  independent  variables x,  ,Xj, .  .  .  ,  .v 
that  minimize  or  maximize  a  given  objective  function  flx^ ,  .v,t ....  xn)  The  objective 
function,  or  performance  function,  as  it  is  frequently  called,  is  repeatedly  evaluated  for 
different  sets  of  independent  variables,  and  each  time  the  function  values  are  compared 
and  the  ‘best'  or  optimizing  variables  are  retained.  The  iterative  manner  in  which  the  inde¬ 
pendent  variables  .tj .  Xj . xn  are  computed  affects  the  speed  and  efficiency  of  the 

direct-search  methods.  In  general,  the  variables  are  changed  ( I )  to  improve  the  value  of 
the  objective  function  and  (2)  to  g>ve  information  about  the  objective  function  which  will 
determine  future  changes  in  the  variables. 

Direct-search  methods  are  analytically  and  computationally  simple  and  h&ye 
proved  successful  in  many  difficult  design  problems.  They  are  also  useful  in  the  solution 
of  t lie  TPBVP  if  it  can  be  formulated  appropriately  as  an  optimization  problem  The 
TPBVP1  is  discussed  in  detail  in  the  following  section  It  is  difficult  to  solve  because  some 
of  the  boundary  conditions  are  specified  at  the  initial  time  while  others  are  specified  at  the 
final  lime.  Thus,  the  solution  requires  that  a  system  of  ordinary  differential  equations  be 
numerically  integrated  until  all  the  boundary  conditions  are  satisfied.  This  problem  there¬ 
fore  requires  that  an  iterative  procedure  be  employed.  As  with  all  iterative  procedures, 
convergence  is  an  important  consideration. 

The  TPBVP  occurs  commonly  in  many  of  the  physical  sciences.  It  occurs  in 
flight  mechanics  for  orbit  determination  as  well  as  in  intercept  and  rendezvous  studies 
It  also  occurs  in  science  involving  the  transport  and  distribution  of  mass  or  energy  Other 
examples  are  the  stress  and  strain  distribution  in  a  given  material  under  a  specified  load 
and  the  flow  of  material  through  a  heat  exchanger  Perhaps  the  greatest  emphasis  on  the 
TPBVP  has  occurred  in  the  study  of  optimal  control  theory.  In  minimizing  a  time  integral 
performance  function  subject  to  the  differential  equations  describing  a  given  system,  a 
TPBVP  must  be  solved.  The  major  methods  of  TPBVP  solution  are  discussed  and  evaluated 
in  a  recent  article  by  Taplev  and  Lewallen*  and  their  results  are  used  in  the  comparison 
with  the  direct-search  methods  presented  in  this  report 

After  a  cursory  discussion  of  the  TPBVP,  the  techniques  that  are  employed  are 
discussed.  These  include  the  method  of  Hooke  and  Jeeves3  and  the  method  of  Flood  and 
Leon4 .  A  method  by  Zangw ill5  ,  which  is  an  improved  version  of  Powell's  method6 .  is 
also  presented.  Though  other  search  schemes  are  available,  these  algorithms  are  well  known 
and  demonstrate  the  principle  of  solving  TPBVP  by  using  direct-search  methods. 

These  techniques  arc  applied  to  an  Earth-Mars  trajectory  problem7  which  results 
from  a  minimum-time  optimal  control  problem.  This  control  problem  has  been  used  as  a 
standard  TPBVP  for  trying  the  various  methods  of  solution.  The  techniques  arc  compared 
for 

1 .  Simplicity  of  formulation 

2.  Convergence  time 

3.  Computational  requirements 

4.  Sensitivity  to  initial  conditions 

Numencal  results  of  applying  direct  search  to  the  trajectory  problem  are  tabulated. 


1  See  RKFt'Rt\CHS 
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TWO-POINT  BOUNDARY  VALUE  PROBLEM 


The  two-point  boundary  value  problem  is  stated  as  follows:  given  the  system  of 
differential  equations 


y('>  =  Ay(')) 


where  y(r)  is  an  rt -dimensional  vector 


y1(0  =  [y’jlO.Vjfr), . . .  (2) 

and  where  superscript  T  denotes  matrix  transpose  and  t  is  'he  independent  variable.  The 
initial  conditions  are  given  by  the  p-dimensional  vector  function 

gT(y(r0),r0)  =  [«,(y00),  r0l,f2(yU0),  f0>,  ,gp(yUp),  <‘0>] 


and  the  final  conditions  by  the  q-dimensional  vector  function 

hT(y(rf),  rf)  =  [h,(y(rf),/f),  h2(y(rf),  ff),...,/iq(y(ff),  rf)] 


with  rQ  specified,  q  =  n  +  1  p.  The  assumption  is  that  1(,  the  final  time,  is  not  specified: 
hence,  it  must  be  determined  along  with  the  n-p  initial  conditions.  If  all  the  initial  condi 
tions  g  and  t(  are  known,  the  problem  is  an  initial-value  problem  and  is  solved  directly  by 
integrating  the  system  represented  by  expression (1).  However,  in  the  TPBVP  in  which  cer¬ 
tain  initial  conditions  have  been  replaced  with  final  conditions,  the  solution  of  the  problem 
is  characterized  by  the  use  of  successive  approximations  to  match  the  given  final  conditions. 

An  example  of  a  simple  nonlinear  two-point  boundary  value  is 

Vj  =  >’2  (0  (>) 

=  -  !>',  (0  l  (6> 


W  = 0  <7> 

y2Or)=B  (bi 

where  r()  and  rf  are  given  The  initial  slope  is  unknown,  and  instead  the  final  condi¬ 
tion  is  specified.  The  existence  and  uniqueness  of  the  solution  are  dependent  on  the 

number  B .  The  methods  of  solution  reviewed  are  essentially  computational  realization-ol 
existence  proofs.  Other  questions  arise:  Does  the  problem  have  a  solution1  One  solution 
only'.'  More  than  one  solution9  Reference  1  contains  a  lucid  discussion  of  the  existence 
and  uniqueness  of  solutions  to  T'"  V?. 

Most  of  the  methods  oi  .olving  the  TPBVP  can  be  classed  in  one  of  two  groups 
'  , )  those  that  solve  a  sequence  of  nor, linear  initial-value  pro'  , is  and  (2)  those  that  solve 

sequence  of  linear  boundary  value  problems. 
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‘SHOOTING’  METHODS  OF 
SOLVING  THE  TPBVP: 
SUCCESSIVE  APPROXI¬ 
MATIONS  OF  INITIAL 
CONDITIONS  LEAD  TO 
A  MATCHING  OF  THE 
GIVEN  FINAL  CONDITIONS 


OTHER  METHODS  SOLVE 
A  SEQUENCE  OF  LINEAR 
BOUNDARY  VALUE 
PROBLEMS  CONVERGING 
TO  A  SOLUTION 


1 .  Methods  of  the  first  class  are  called  ‘shooting’  methods.  These  methods  select 
trial  initial  conditions,  integrate  the  system  of  differential  equations,  and  then  check  to  see 
whether  the  final  boundary  conditions  are  satisfied.  The  errors  at  the  final  boundary  are 
used  to  correct  the  initial  values  and  the  process  is  repeated.  In  the  context  of  this  report 
the  shooting  method  is  summarized  as  follows; 

If  the  unknown  initial  conditions  are  identified  as  scalars  Xj,  the  dependence 
of  the  solution  upon  the  initial  conditions  can  be  emphasized  by  writing >’(f,x). 

Therefore,  the  final  boundary  conditions  (4)  will  depend  on  xT  =  (Xj ,  .  .  .  , x  ) 
as  shown  by 

h(x,  rf)  ~  0  (9) 

The  shooting  method  therefore  is  a  method  of  finding  a  x°  solution  to  (9);  that  is, 

h(xu,/f)=0  (10) 

If  the  final  time  rf  i-  unknown,  an  additional  scalar  x  4  ,  =  tf  is  defined  and  equation  (9) 
becomes 


hlx)  =  0  (II) 

with  x1  =  (Xj ,  .v, . x  +( ).  Any  standard  procedure  for  finding  roots  of  a  vector 

function  can  be  used.  Integration  of  the  system  of  differential  equations  is  of  course 
required  io  evaluate  h!x). 

2.  The  other  class  of  methods  for  solving  TPBVP’s  involves  solving  a  sequence  of 
linear  boundary  value  problems  such  thal  the  sequence  converges  to  the  solution  of  the 
original  problem.  Among  these  methods  are  quasilinearization  and  the  method  of  perturba¬ 
tions.  The  methods  are  useful,  since  solving  linear  TPBVP's  is  considerably  simpler,  though 
formulating  and  computing  with  these  methods  are  more  complex  operations.  The 
method  of  perturbations  is  presented  to  illustrate  this  class  of  methods. 

The  method  of  perturbations  develops  a  Taylor’s  series  expansion  of  the  original 
problem  and  solves  the  linear  problem  that  results.  The  >r(h  approximation  is  g;,cn  by 


y*+t  =  f.i%9  iy„*rv„.]  +°2(y„try»> 


(121 


If  the  notation  fiy  =  v  ,  -  v  is  used  and  onlv  the  linear  terms  in  the  expansion  are  re- 
1  ■  n*  \  •  r  ■  r 

tamed,  the  linear  perturbation  is  seen  to  be 


6y 


6y  =  A5y 


where  the  matrix  of  the  partial  derivatives 


(13) 


A  = 


is  evaluated  on  the  itth  trajectory  (note  that  integrating  the  nonlinear  equations  is  still 
required).  The  final  boundary  conditions  are  also  expanded  about  the  terminal  conditions 
of  the  ii'1’  trajectory.  The  result  can  be  expressed  as 


<U) 
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When  ihe  effect  of  a  variarion  in  the  terminal  time  is  considered 

dy/ =  *y/+y/d,f  (15) 

the  final  condition  (14)  becomes 

dh  =  jjpjfiy,  +  hdtf  (16) 

Solving  equation  (13)  puts  the  terminal  conditions  in  terms  of  the  initial  condition  varia¬ 
tions  (6yj  -  irSy^)  where  tr  is  the  characteristic  solution,  and  the  above  becomes 

dh  =  [i7>6yo+hd'*  (,7) 

Furthermore,  since  the  variations  in  y0  must  satisfy  g(y(f0),  /0)  =  0,  follows  that 

d,‘[|]  6y»'°  <l8> 

*0 

Hence,  once  dh  is  determined,  equations  ( 17)  and  ( i8)  are  solved  for  unknown  6y0  and 
dq  .  The  procedure  is  then  repeated  until  dh  is  within  some  specified  convergence 
criterion  about  0. 

The  methods  in  this  category  are  many,  with  various  differing  analytical  and  com¬ 
putational  properties  These  methods  are  discussed  briefly  in  order  that  comparisons  can 
oe  made.  The  method  of  perturbation  is  recognized  as  a  simple  method  to  formulate.  It 
is  computationally  simple,  requiringonly  the  integration  of  the  original  differential  equa¬ 
tions  to  evaluate  the  linearized  equations 

The  quasilinearization  method  is  given  by  the  equation 


where  n  is  the  iteration  number  here.  This  equation  is  solved  repeatedly.  Past  informa¬ 
tion  is  used,  and  both  y  and  y  (for  each  n)  from  ihe  previous traieciorv  are  stored  The 

n  hj  ^ 

storage  problem  is  severe.  The  results  generally  are  good.  (See  Tapley  and  Lewallen*  and 

Kenneth  and  McGill  K' 

In  this  study  the  above  methods  were  not  actually  tried.  Results  given  by  Tapley 
and  Lewallen  are  used  in  the  comparison  with  the  direct-search  approach  to  TPBVP 
that  follows. 


DIRECT-SEARCH  TECHNIQUES 


Tlte  theory  of  direct  search  is  employed  to  find  a  solution  to  the  unconstrained 
optimization  problem 


min/(x)  (20) 

x 


where  x  e  E"  is  an  rr-dimensional  Euclidean  vector  space  with 


DIRECT  SEARCH 
METHODS  LEND 
THEMSELVES  TO 
OPTIMIZATION 
PROBLEMS 


HOOKE  AND  JEEVES 
USED  AN  'EXPLORATORY 
MOVE' OF  SMALL. 

FIXED  SIZE.  FOLLOWED 
BY  A  'PATTERN  MOVE' 


but  without  using  the  usual  necessary  conditions  lor  an  optimum;  for  example,  V/lxn)=  0, 
or  higher-order  derrvative  information.  The  vector  variable  x  ts  changed  in  a  systematic 
manner  and  the  scalar  objective  function  evaluated  after  each  change.  The  values  of  the 
objective  function  are  compared  and  changes  in  the  variables  are  continued  until  a  value 
of  x  is  found  which  minimizes /fx).  The  theory  in  changing  (he  variables  affects  the 
analytical  and  computational  properties  of  these  methods,  and  these  techniques  are  dis¬ 
cussed  in  the  ensuing  paragraphs.  In  general,  these  methods  are  easily  understood  and 
easily  applied.  Since  they  do  not  require  the  computation  of  the  gradient,  they  lend 
themselves  to  optimization  problems  in  which  the  gradient  is  extremely  complex  and 
difficult  to  compute  or  the  objective  function  is  discontinuous. 

The  difficulty  in  solving  the  optimization  problem  (20)  is  caused  by  the  fact  that 
the  surface  represented  by 

z  =  /lx)  (22) 

is  often  not  of  a  convenient  'shape.'  It  would  be  convenient  if /fx)  had  a  bowl  shape, 
with  its  minimum  clearly  defined.  Then  a  search  through  various  x  would  yield  the 
optimum.  However,  the  surfaces  often  contain  troughs,  and  direct-search  methods  often 
move  about  inefficiently  looking  for  the  optimum.  The  philosophy  of  search  must 
include  a  technique  to  seek  out  a  trough  or  ridge  and  follow  it  to  an  optimum.  Each  of 
the  search  techniques  that  follows  has  this  property  to  some  degree. 

The  method  of  Hooke  and  Jeeves  was  one  of  the  earliest  and  most  accurate  of 
the  search  methods.  It  possesses  the  ability  to  follow  a  ridge,  a  properly  that  accounts 
for  its  success.  Further,  Hooke  and  Jeeves  introduced  certain  terms  which  have  become 
standard  terminology  in  search  methods.  The  ‘exploratory  move’  is  a  systematic  search 
of  the  coordinate  directions  to  determine  an  improved  point.  The  method  of  Hooke  and 
Jeeves  (see  illustration)  performs  an  exploratory  move  with  a  small,  fixed  step  size. 

Mathematically,  the  exploratory  search  evaluates 

JixrX2, - x±A,. ..,*.)  t=l,2,.  ,ri  (23) 

and  compares  it  with  the  best  optimum  found  on  the  previous  coordinate  search.  This 
new  point  found  as  a  result  of  an  exploratory  move  is  called  a  base  point .' 

The  'pattern  move’  always  follows  a  successful  exploratory  move.  The  direction 
for  a  pattern  move  is  determined  by  the  previous  base  point  xB  and  the  best  point  of  a 
successful  exploratory  move  -  The  dnection  is  given  by 

d  =  x(  -  xB  (24) 

The  point  Xj. .  the  result  of  the  last  exploratory  move,  becomes  the  new  base  point  and  a 
new  exploratory  move  is  made  from  x!  +  d.  As  long  as  the  exploratory  move  yields  an 
improved  point,  the  distance  x^  -  xB  increases  and  the  method  accelerates  and  simultane¬ 
ously  follows  a  ridge.  When  an  exploratory  move  fails,  the  step  size  in  the  search  scheme 
is  reduced  by  a  reduction  factor  (<  1 .0)  and  the  method  continues. 

This  method  does  not  determine  the  direction  of  steepest  descent  but  settles  for 
any  improvement.  However,  its  choice  of  direction  yields  an  efficient  and  accurate 
search  method.  It  is  implemented  in  a  computer  routine  called  NELC  DIRECT  1 

The  method  of  Hooke  and  Jeeves  (NELC  DIRECT)  has  no  formal  convergence  i 

theorems  to  support  its  success,  but  it  works  and  works  well.  1 
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FLOOD  AND  LEON  USED 
EXPLORATORY  MOVES 
TO  FIND  A  MINIMUM 
ALONG  COORDINATE 
DIRECTIONS,  THEN  A 
PATTERN  MOVE 


The  method  of  Flood  and  Leon  (UN1VAR)  also  has  no  formal  theory  but 
appears  intuitively  to  be  more  plausible.  It  makes  an  exploratory  move  that  searches 
along  each  coordinate  direction  until  a  minimum  is  found.  Then  it  makes  a  pattern  move 
using  the  original  base  point  and  the  final  improved  point  of  the  exploratory  move  to 
determine  the  new  direction.  The  method  then  minimizes  along  that  direction  In  mini¬ 
mizing  along  a  given  direction,  the  method  accelerates  by  taking  increasingly  larger  steps. 
The  method  is  briefly  described  by  the  following  steps: 

1 .  The  exploratory  move  is  made  by  incrementing  an  element  of  the  vector  x,  such 

that 

J\x , ,  x , . X  +  pkA _ X  )<fix..x, . x.  +  p.k'l± . x 

J  I  2  ■  i  *»  ’  '  it  •'  1  2  i  i  n  1  ' 


is  a  minimum,  p.  is  a  number  greater  than  I ,  and  by  raising  it  to  successive  powers, 
is  successively  increased,  resulting  in  an  acceleration  along  that  direction.  Successive 
moves  are  continued  until  a  failure  occurs  (functional  value  gr eater  than  the  previous 

value),  then  the  procedure  sets  x  =  ar. , .  .  .  ,.v  *  pk  A . x  and  reinitializes  k  to  zero. 

Two  successive  iterations  of  this  procedure  are  made  to  find  a  minimum  along  the 
coordinate  direction.  If  a  success  is  not  found  with  k  =  0,  the  direction  ol  search  along  the 
coordina.e  directions  is  reversed,  that  is,  xv  -  pk  A.  If  this  fails,  that  point  is  assumed  to 
be  the  minimum. 


2.  A  pattern  move  is  made  if  the  functional  value  at  the  end  of  the  exploratory 
move  is  better  than  the  point  at  the  beginning  of  the  exploratory  move.  The  direction  of 
the  move  is  determined  by  the  difference  of  the  last  point  x(  and  the  beginning  point  xB 
of  the  exploratory  move.  It  is  given  by 

x  -  x,.  +  pk  (x,  -  xHl  (26) 

and  the  accelerating  factor  p  is  raised  to  the  A11’  power  I  k=  1 .  2 ....  I  so  as  to  accelerate 
along  that  direction  (A  repi events  the  number  of  successes  in  a  given  direction).  After 
failure  occurs,  the  procedure  outlined  in  step  I  is  repealed  resetting  k  to  zero.  The  nega¬ 
tive  direction  is  searched  if  failure  occurs  at  A  =  0. 

The  method  of  Zangwill  is  a  search  method  with  proved  quadratic  convergence 
It  is  described  as  follows: 

Let  c  ,r  =  1 , 2,  ...  ,ii  be  the  unit  coordinate  directions.  The  coordinate  direc¬ 
tin'’'  are  used  after  one  iteration  of  the  following  procedure:  ( I )  Beginning  with  an  initial 
p  ii. i  and  n  directions^  .  the  method  minimizes  /(x  ,  +  a i .)  selling  x  =  x  ,  +  a  £ 

r  0  i-l  i  ;  -  i  i-l  i”: 

foi  i  =  1 ,  2,  .  . . ,  n.  (2)  A  direction  +  j  is  computed . 

X  -  xn 

=  nx'  -  x0i:  ,27) 

where  ||x(j  -  x0 1|  is  the  Euclidean  norm  (square  root  of  the  sum  of  the  squares  of  the 
elements)  and  x  ..  =  x  +  a-  ,£  „ , ,  where  a  , ,  minimizes /(\  +  a£  ,  ). 

’  1  II  11+  1^11+  I  il+  1  n  ^ji+ 1 

The  directions  are  updated  as  follows: 

4,=4,+  1  (=l,2.--.,u  (28) 

is  removed  from  the  basis.  (3)  A  direction  c(  is  introduced  and  the  point 
x  .-,=x  ,,+q  ..c  is  determined  so  that  rtx  , .  +  a.  ,c  )  with  respect  to  a  ,  is  a 
minimum.  It  +  +=  0  the  index  r  (f  =  1 ,  2, .  .  -  , is  updated  and  the  routine  goes  to 

part  ( I )  If  a  2  =  0  for  /  =  1 . n  the  solution  has  been  found.  Note  that  r  is  cycled 
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POSING  THE  TPBVP  AS 
AN  UNCONSTRAINED 
OPTIMIZATION  PROBLEM 
NECESSITATES  FORMU¬ 
LATING  AN  OBJECTIVE 
(OR  WEIGHTING)  FUNCTION 


through  the  coordinate  directions;  it  is  incremented  through  n  and  then  reset  to  I .  This 
prevents  the  directions  4^,  k  =  1 , 2, .  .  .  ,  n  from  becoming  linearly  dependent.  This 
routine  is  shown  to  yield  convergence  for  a  quadratic  in  a  finite  number  of  iterations. 


FORMULATING  A  MATHEMATICAL  PROGRAMMING 
PROBLEM 


The  two-point  boundary  value  problem  must  necessarily  be  posed  as  an  uncon¬ 
strained  optimization  problem  before  the  direct-search  methods  can  be  applied.  The 
formulation  of  the  unconstrained  optimization  problem  can  be  obtained  by  an  extension 
of  the  shooting  method  (see  TWO-POINT  BOUNDARY  VALUE  PROBLEM)  The 
terminal  or  final  conditions  (4)  for  the  system  of  differential  equations  are  specified,  and 
for  any  given  set  of  initial  conditions,  the  terminal  conditions  are  computed  and  com¬ 
pared  with  the  specified  terminal  conditions.  The  goal  is  to  make  the  difference  between 
the  specified  and  computed  final  boundary  conditions  as  small  as  possible;  therefore,  an 
objective  (or  ‘weighting,’  or  ‘penalty’)  function  of  the  differences  is  minimized  by  choos¬ 
ing  the  unknown  initial  conditions.  Because  the  specified  final  boundary  conditions  are 
zero  in  the  TPBVP,  the  objective  function  must  have  the  following  properties: 

F{ h)  =  0  when  h(y(rf  >,  rf)  =  0 
and 

F(h)  >  0  when  lh(y(/f). /()|  ¥■  0  (29) 

One  posable  form  of  F ,  the  objective  function,  is  the  sum  of  the  absolute  values  of 
h'iy'jf ),  /f);  namely 


Q 

F(h)  =  |F(x)|  =  F(x)  (30) 

/  1 

where  x  replaces  those  initial  conditions  that  are  unknown  and  xn  =  r( .  Hence,  the 
weighting  function  F  is  a  function  of  the  initial  conditions  and  the  final  time  t(,  and  the 
problem  is  now  in  the  form  of  an  unconstrained  minimization  problem 

niinF(h(x))  =  min  F(x)  (31) 

x  x 

The  various  direct-search  algorithms  are  now  applicable.  It  is  understood,  however,  that 
to  evaluate  F(x)  we  must  integrate  a  system  of  nonlinear  differential  equations  with  the 
unknown  initial  conditions  and  final  time  r,  specified  as  x  by  the  algorithm. 

Another  form  of  weighting  function  which  is  used  often  in  practice  and  which  is 
employed  in  the  trial  problem  is 

F(x)  =  h(x)T  h(x)  (32) 

This  is  the  sum  of  the  squares  of  the  errors  in  the  final  conditions,  and  the  malhemalicai 
programming  problem  takes  the  form  of  least  squares,  that  is,  x  is  chosen  so  as  to  mini¬ 
mize  the  sum  of  the  squares  of  the  error. 

Other  independent  constraints  may  be  imposed  in  the  formulation  of  the 
TPBVP,  but  these  have  not  been  considered  here.  There  are  certain  constraints,  such  as  on 
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the  variable  t(  (for  example,  tf  >  0),  which  can  be  included  as  a  penalty.  Any  time  t(  vio¬ 
lates  the  constraint,  F  is  equated  to  a  very  large  number  ( 1 06 ).  In  general,  the  addition 
of  other  constraints  will  result  in  a  far  more  complex  problem. 

The  formulation  of  the  TPBVP  as  a  mathematical  programming  problem  is  easily 
effected.  It  remains  to  discuss  how  successful  direct  search  is  in  its  solution. 


RESULTS  AND  CONCLUSIONS 

The  technique  of  using  direct  search  was  applied  to  a  two-point  boundary  value 
problem  that  resulted  from  the  calculation  of  a  minimum-lime  Earth-Mars  trajectory. 

This  particular  problem  has  been  extensively  used  by  many  investigators7  7 ’*  in  studying 
methods  of  solving  the  TPBVP,  and  some  valid  comparisons  can  be  made  between  direct- 
search  methods  and  other  recognized  methods.  The  trajectory  problem  and  a  list  of  the 
constants  are  given  in  appendix  A.  The  results  of  computer  runs  on  the  CDC  1604  are 
tabulated  on  the  following  page. 

The  general  characteristics  of  TPBVP  solution  methods  considered  are:  (1)  sim¬ 
plicity  of  formulation  and  implementation,  (2)  computer  storage  requirements,  (3)  con¬ 
vergence  sensitivity,  and  (4)  convergence  time.  Of  these  characteristics,  the  last  is  the 
most  difficult  to  compare  realistically.  Different  numerical  integration  schemes  as  well  as 
different  computers  may  significantly  affect  convergence  time,  and  any  comparisons  are 
dependent  on  these  variables.  In  contrast,  the  other  characteristics  are  independent  of  the 
numerical  integration  scheme  and  computer. 

1.  Simplicity  of  Formulation 

The  direct-search  method  applied  to  the  TPBVP  is  by  far  the  simplest  to  formu¬ 
late  and  implement.  A  penally,  or  weighting,  function  of  the  terminal  boundary  condi¬ 
tions  is  chosen.  The  calculation  of  the  penally  function  requires  only  the  integration  of 
the  differential  equations  and  is  therefore  simpler  than  any  method  involving  the  calcula¬ 
tion  and  integration  of  additional  perturbation  equations.  The  initial  conditions  are 
specified  directly  by  the  search  method  that  is  employed. 

2.  Computer  Storage  Requirements 

Computer  storage  requirements  in  applying  the  direct-search  methods  to  TPBVP’s 
are  minimal.  Since  only  forward  integrations  of  the  differential  equations  are  necessary , 
no  intermediate  storage  is  required.  Other  methods  for  TPBVP  solutions  are  more  com¬ 
plex,  since  the  intermediate  values  of  differential  equation  variables  must  be  stored  for 
future  integration.  Such  storage  poses  a  severe  problem  in  solving  large  TPBVP’s. 

3  Convergence  Sensitivity 

Methods  of  direct  search  are  highly  insensitive  to  initial  conditions.  Convergence 
sensitivity  of  the  direct  methods  was  tested  with  +20-percent  change  in  the  initial  guess 
of  final  time  .  All  methods  converged.  The  method  of  UNIVAR  demonstrated  the 
capability  of  searching  a  large  parameter  space  in  a  short  time.  NELC  DIRECT  converged 
for  large  changes  m  initial  conditions  -  however,  w  ith  a  larger  computer  run  time.  The 
insensitivity  to  initial  conditions  uf  these  methods  is  a  decided  advantage  of  this  approach 
to  solving  TPBVP’s. 
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VARIOUS  SEARCH  ALGORITHMS  (C 


4.  Convergence  Time 

Convergence  times  of  the  direct  methods  could  not  be  realistically  compared. 

For  different  step  sizes  and  different  search  methods,  different  convergence  times  result. 
The  tendency  in  the  computational  work  was  to  choose  the  step  size  5  too  small  in 
UNIVAR  and  NELC  DIRECT.  Often  the  larger  step  size  results  in  quicker  convergence 
but  less  accuracy.  However,  there  are  results  which  stand  out  for  their  accuracy.  There 
are  other  factors  in  the  UNIVAR  search  scheme,  such  as  the  search  acceleration  factor  p, 
which  further  compound  the  difficulty  of  convergence  evaluation.  The  acceleration  factor 
was  chosen  at  1 .618  and  never  changed  in  the  computation.  Zangwill’s  method  gave  the 
greatest  accuracy  with  the  shortest  run  time  (10  minutes’  maximum  run  time). 

In  general,  the  search  schemes  performed  remarkably  weli.  Simplicity  stands  out 
as  (heir  strongest  computation  characteristic.  The  results  were  satisfactory  with  Zangwiil’s 
search  method,  which  has  quadratic  conveigence  properties.  However,  the  methods  of 
UNIVAR  and  DIRECT  exhibited  dependence  on  the  choice  of  search  step  size  (delta  or 
deltai)  and  the  rho  accelerating  factor.  As  the  theory  of  search  methods  evolves,  additional 
algorithms  will  be  developed  which  will  add  to  the  attractiveness  of  the  direct-search 
approach  for  TPBVP’s.  New  and  faster  computers  will  further  contribute  to  the  utilization 
of  these  methods  for  solution  of  the  TPBVP.  However,  the  direct-search  approach  always 
required  ionger  computer  time  to  reach  a  solution  than  the  indirect  methods  given  in 
references  2,  8,  and  7. 


REFERENCES 


1 .  Bailey,  P.  B.  and  others.  Nonlinear  Two  Point  Boundary  Value  Problems, 

Academic  Press,  1968 

2.  Tapley,  B,  D.  and  Lewallen,  J.  M.,  “Comparison  of  Several  Numerical  Optimization 

Methods,”  Journal  of  Optimization  Theory  and  Application^.  l,p.  1-32, 
July  1967 

3.  Hooke,  R,  and  Jeeves,  T.  A.,  “  ‘Direct  Search’  Solution  of  Numerical  and  Statistical 

Problems,"  Association  for  Computing  Machinery.  Journal,  v.  8,  p.  212-229, 
April  1961 

4.  California.  University.  Space  Sciences  Laboratory  Internal  Working  Paper  20, 

Comparison  Among  Tight  Known  Optimization  Procedures,  by  A  Leon, 
August  1964 

5.  Zangwill,  W.  1.,  “Minimizing  a  Function  Without  Calculating  Derivatives,” 

Computer  Journal,  v.  10,  p.  293-296,  1967-1968 

6.  Powell,  M.  J.  D.,  “An  Efficient  Method  For  Finding  the  Minimum  of  a  Function  of 

Several  Variables  Without  Calculating  Derivatives,"  Computer  Journal,  v.  7, 
p.  155-162,  1964-1965 

7.  Kelley,  H.  J.,  “Method  of  Gradients,”  Chapter  6  in  Optimization  Techniques,  by 

G.  Leitmann,  Academic  Press,  1962 

8.  Kenneth,  P.  and  McGill,  R.,  “Two-Point  Boundary-Value-Problem  Techniques,” 

Chapter  2  in  Advances  in  Control  Systems,  v.  3,  Academic  Press,  1966 


APPENDIX  A:  THE  OPTIMAL  CONTROL  PROBLEM 


The  problem  is  to  determine  the  minimum-time  trajectory  from  Earth  to  Mars. 
The  optimal  control  problem  is  to  minimize 


subject  to 


(*2)"  K.  r 
x ,  - - +  —  sin  8 

1  *3  U3)2  m 

•  *2  T  n 

X-,  =  -  +  —  cos  0 


**3  ”  A1 

with  respect  to  x,  X,  rf,  and  0  (these  equations  are  derived  in  appendix  B),  with 

m  -  m0  +  ct  (t 

x,  (rQ),x,  (r  l.Xj  (/q),  and  x,  (/f  l.x^  (/f),Xj  (ff)  given  (K,  T,  m0,  and  c  ate  known 
constants)'. 

The  problem  can  be  solved  by  using  the  maximum  principle.  A  generalized 
Hamiltonian  is  formed 


H  =  -1  +X 


(*2)2  K.  T  1  .  f  *t  *2  T  1 

L  *3  (x,)2  m  J  1 L  *3  m  J 


+  X,  X.  (A-6) 


The  necessary  conditions  are 


z  3W  a2  *2  , 

x‘ =  ■  5r; =  x3  ■ 

(A-7) 

..  3//  2jc2  X,  .x,  \2 

2  3x2  xi  x3 

( A -8) 

3 H  M*2)2  ,  K  X2  r2  *1 

dxi  (X3)2  (x3)3  (x3)2 

(A-9) 

(-v2>2  k  r 

x.  -  -  -+  sin  £ 

’  v3  (x3)2  m 

(A- 10) 

xlx2  T  „ 

x+  =  -  +  cosp 

-  m 

( A- 1  1 ) 

x3  -  X  j 

(A-12) 

bH  T  ,  .  T  .  .  r 

TT  =  A »  —  cos  p  -  X,  —  sin  =  0 

3/3  ■  /n  2  m 

(A- 13) 

x(  (/y),x7  (/  >,x3  |f  ),X|  (f,),x2  (f,).x.  (rf)given  with  transversality  condition 
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+  —  sin  (3 )  +  X, 


.+_c°sij  +  X3.v, 


=  0 


(Further  discussion  of  this  problem  is  found  in  reference  1 .) 
It  is  possible  to  solve  explicitly  for  (3 

X, 

tan  ,3  =  t— 

A2 


(A- 14) 


(A- 15) 


or 


sin  (3  = 


cos  d  = 


n/(X,)2-MA2)2 

)2  +  (X2>2 


The  transversality  conditions  (A-14)  can  be  modified  to 

j-  (X,  (r,))2  +  (X2  (rf))2 
V(X,  (ff))2  +(X2  <?f)>2 

-  I  +  £  x/(X,)2  +(X2)2  =  0 


x/(Xj  )2  +  (Xj  )2  =y 


( A-16) 

(A  17) 


(X,)2  +  (X2 12 


m  * 

r 


( A-18) 


where 


(fr)  =  0, 


<*2  (ff)r 


— S_1 

Uj  (t,-))2  J 


In  minimum  time  problems,  one  multiplier  can  be  specified  and  hence  the  number  of 
unknowns  reduced.  This  occurs  when  the  state  equation  depends  only  on  the  ratio 
Xj/X2,  and  it  is  possible  to  replace  the  Final  condition  (A-18)by  a  normalized  variable 
Xj  (rQ)  =  1  0.  The  transversality  condition  is  eliminated  and  the  number  of  unknown 
initial  conditions  is  reduced  by  one. 

The  TPBVP  is  now  given  by: 


<  K  ^  T  X) 

V3  (,r3)2  \/ (X, )-  +7xT7 


( A-19) 


*1  *;  T_  X-, _ 

xi  m  VfX,  )2  +  (X2r 

*3  =  -V, 

X2.x2 


(A-20) 


( A- 2 1 ) 


(A  22) 


Ik  »la."-k  Ha.ifc  -„.rfi.-4lMiri»lli4i^ill-iMMi.<.-,-4ll|J>illrt*«|iliii  MlM.  •  >.*  ul  I.IrtJ- ,»  i,  l.«  -  "  >44«:ilil>j  Hjilai  "  .Ullftirit-*  II  iwlui*  aJIliatiM’llU'iittltWia1  •  .  i H>  ■4»">iH<*luilC4r..  a*«ll*W.W*.i*.ii.-  odbMLfclhUHM,  jiMhiAULl 


iJf2  *1  X\  *2 

X,  =  -  + 

2  *3  *3 

(A-23) 

X,  (*2)2  2K  X2JC2Jfl 

(x3)2  (,r3)3  U3)2 

(A-24) 

‘-C 

o' 

1! 

O 

o 

(A-25) 

x2  (,o)  =  10 

(A-26) 

x3  (/Q)  =  1.0 

(A-27) 

X,  (/„»  =  1.0 

(A-28) 

x,  <ff)  =  0.0 

(A-29J 

.v,  <rf)  =  0.8098 

( A-30) 

,v3  (/,-)  =  1.525 

( A-3 1 ) 

in  the  FORTRAN  routine  with 

|.V|  . x,,x3.  Xj ,  X3 .  X, ] 

( A-32) 

The  computer  routine  is  self-contained:  comment  cards  indicate  the  necessary  data  to  be 
specified  (appendix  C). 


APPENDIX  B:  DYNAMICS  FOR  A  TRAJECTORY 
OPTIMIZATION  PROBLEM 


The  problem  to  be  solved  involves  the  control  of  a  constant-thrust  rocket  en 
route  front  Earth  to  Mars.  The  angle  of  thrust  with  respect  to  a  local  horizon  (see  page  20) 
is  the  control  variable.  The  problem  is  assumed  coplanar  (that  is,  Earth  and  Mars  are 
assumed  to  be  points  in  a  plane)  with  Earth  and  Mars  exhibiting  negligible  gravitational 
force.  The  derivation  of  the  dynamics  involves  Newton’s  law  (F=  ma).  In  a  polar 
coordinate  system,  Newton’s  law  takes  the  following  form: 


/  d‘0  ,  dr  d0  \  v , 

Ir —  +  0 - I*!/, 

\  dr  drdf/ 


(B-l) 


(B-2) 


where 


—  =  u  =  radial  velocity 
dr 

d0  .  ,  , 

r —  -  v  -  tangential  velocity 
dr 


.  .  d0  dz6  d0  d~0 

v  -  r—  +  r -  =  u  —  +  r - 

d'  dr2  d'  dr2 


/  are  the  forces 

r  is  the  radial  position  of  the  vehicle 
0  is  its  angular  position 

The  first  equation  (B-l)  of  motion  is  rewritten 


/•  d0  d0  \ 

jlu-  v~\  -  -H;-7'sinJJ 


where 


and 


since 


W  =  weight  of  the  rocket 


T  =  thrust  of  the  rocket 


d0  _  v 
dr  r 


(B-3) 


(B-4) 


(B-S) 


Equation  (B*4)  becomes 


■c 


-  b'  +  T  sin  P 


u 


r 


H-'  T  .  „ 

—  +  —  sm  (3 

m  m 


The  second  equation  of  motion  becomes 


(B-6) 


or 


*  u  v  ~  a 

m  v  +  rri  —  '  T  cos  p 
r 


<B-7) 


which  is  written 

•  uv  .  V’ 

p  =  —  +  —  cos 
r  m 

One  further  equation  is  required  -  the  equation  describing  change  of  mass  with 


time.  This  equation  is  given  by 


m  =  mQ  +  cf 

(B-8) 

which  assumes  a  constant  change  of  mass. 

These  equations  arc  a  simplification  of  the  trajectory  problem  but 

approximate  model  of  the  system. 

The  state  variable  model  is  obtained  by  defining 

represent  an 

*1 

=  u  =  —  =  radial  velocity 
dr 

(B-9) 

x2 

=  v  =  tangential  velocity 

(B-10) 

xi = 

r  =  radial  position  of  vehicle 

(B-l  1) 

The  equations  are 

*1 

K  J  .  „ 

=  — : - +  —  sin  p 

XJ  (Xj)'  m 

(B-l  2) 

x  ,  * - +  —  cos  P 

x}  m 

(B-l  3) 

*3  =  *1 

(B-l  4) 

where 

vv  _  K 

(x,>2 

and  K  is  the  gravitational  constant  of  the  sun. 


19 


min»w  !ipi  *  mppn  nomiM;^ . hr  i|v  w  i»i  '■  ■"•  "ffW"  l-T'-q  ■  'Wll1f>fM,l»'''inr  Iff  WW  ll*fJ|R,'lWwl 


1 


A 


APPENDIX  C:  COMPUTER  FLOW  CHARTS  AND  INSTRUCTIONS 


1969  045 

PROGRAM  TPBVPSCH 

C  THIS  PROGRAM.  WILL  SOLVE  A  TWO-POINT  BOUNDARY  VALUE  PROBLEM 
C  USING  A  DIRECT  SEARCH  SCHEME  (UNIVAR)  AND 
C  THE  INTEGRATION  SCHEME  IS  THE  FOURTH  ORDER  RUNGE-KUTTA 
C  AS  MODIFIED  BY  GILL 
C 

C  THE  XI  ARE  THE  INITIAL  CONDITIONS  THAT  ARE  UNKNOWN  AND  WHICH 
C  ARE  SPECIFIED  INITIALLY  IN  THE  PROGRAM  DRIVER 

C 

C  K  REPRESENTS  THE  NUMBER  OF  UNKNOWN  VARIABLES. 

C  N  =  THE  NUMBER  OF  DIFFERENTIAL  EOUATIONS  TO  BE  INTEGRATED 
C 

C  THE  DELTAI  REPRESENTS  A  STEP  CHANGE  IN  THE  INITIAL  CONDITIONS 
C  AND  WHICH  IS  USED  BY  THE  SEARCH  ROUTINE. 

C  THE  DELTAI  SHOULD  BE  LARGE  INITIALLY.  THEN  IF  RESULTS  ARE  MOT 

C  SAT  IS FACTORY, DELTAI  CAN  BE  REDUCED. 

C 

C  THE  SUBROUTINE  FN  IS  USED  FOP.  INTEGRATION  AND  ALSO  TO  COMPUTE 
C  A  WEIGHING  FACTION  OF  THE  FINAL  BOUNDARY  CONDITIONS 
C 

C  A  WEIGHING  FUNCTION  OF  FINAL  CONDITIONS  MUST  BE  FORMULATED 
C  THE  WEIGHING  FUNCTION  IS  USED  BY  THE  DIRECT  SEARCH  SCHEFE 

C  TO  DETERMINE  THAT  THE  FINAL  CONDITIONS  ARE  MET 

C  THE  FINAL  CONDITIONS  ARE  WRITTEN  AS  EOUALITIES  P(Y)=0.0 

C  THE  WEIGHING  FUNCTION  IS  THE  SUM  OF  THE  SOUARES 

C  G=P(YCl)>PCYO)>+PCYC2)>P(Y(2))  4  ETC. 

c 

C  IF  THE  FINAL  TIME  IS  UNKNOWN  ,  AN  ADDITION  VARIABLE  XICK  )  IS 
C  USED  TO  REPRESENT  THE  FINAL  TIME  AND 

C  INSERTED  IN  THE  INTEGRATION  STOPPING  CRITERION .  OTHERWISE 

C  THE  TRUE  FINAL  VALUE  OF  THE  INDEPENDENT  VARIABLE  IS  USED 
C 

C  PROBLEMS  CAN  BE  ENCOUNTERED  WITH  STEP-SIZE  DELTAI  AND  THE 
C  INTEGRATION  STEP-SIZE  H.  THE  LOCAL  VS  GLOBAL  PROBLEM  IS  PRESENT 
C  OCCURRED  AND  IT  IS  NECESSARY  TO  RESTART  THE  ROUTINE 

C  WITH  NEW  INITIAL  CONDITIONS  TO  GUARANTEE  A  GLOBAL  OPTIMUM  SOLUTION. 

C 

C  ADVANTAGE  OF  THIS  ROUTINE  IS  THAT  IT  ONLY  REQUIRES  INTEGRATION 

C  OF  THE  DIFFERENTIAL  EQUATION  .  NO  PERTURBATION  EOUATIONS 

C  ARE  REQUIRED. 

C 

DIMENSION  XIC4),X(4),DPC4),DELC4),XRC4) 

EXTERNAL  F 
C 


XIC1)=0.7 

XI(2)=0.3 

XIC3)=3.03 
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C 

C  K=  NO.  OF  UNKNOWNS 


C  INITIAL  DELTA 

DELTAI=0.04 

RH0=1.618 

CALL  UNIVAR(K,XI,RHO,DELTAI,X,DP,DEL,XR,FR) 
END 


1969  043 

SUBROUTINE  FNCG,K,XI) 

DIMENSION  XICK) 

C 

C  THESE  QUANTITIES  ARE  NEEDED  FOR  INTEGRATION.  THE  ARRAYS  HAVE  THE 
C  SAME  ORDER  AS  THE  DIFFERENTIAL  EQUATIONS 
DIMENSION  Y(6),FC6) 

TYPE  REAL  MAGSO,MAG,MOM 
C 

EXTERNAL  RKLDEQ 
N=6 
C 

C  THE  KNOWN  INITIAL  CONDITIONS  ARE  SPECIFIED  HERE 
Y(1)=0.0 
Y(2)=1.0 
YC3)=1.0 
YC4)=XIC1) 

YC5)=XIC29 

Y(6)=1.C 

C 

C  THE  INITIAL  VALUE  OF  X,,  THE  INDEPENDENT  VARIABLE  AND  THE  INTEGRATION 
C  STEP  SIZE  H  ARE  SPECIFIED  .  THE  RUNGE-KUTTA  INTEGRATION  SCHEME 
C  IS  SENSITIVE  TO  STEP  SIZE  AND  H  MUST  BE  CHOSEN  CAREFULLY 
H=0 .030$  X-0.0 

c 

2  NT=0 
C 

C  THE  RIGHT-HAND  SIDE  OF  THE  DIFFERENTIAL  EQUATIONS 
C  DY/DX  =  FCY,X) 

C  IS  EVALUATED  NEXT 
C 

6  MAGSQ=Y(4)*Y(4)+Y(5)’!!Y(5) 

QM=1.0 

MAG=SORTFCMAGSO) 

MOM=1.0-0.0749*X 
SB=  Y(4)/MAG 
CB=  Y(5)/MAG 

FC1)=YC2),<‘YC2)/YC3)-GM/(yC3)^YC3)>0. 140500  *SB/M0M 
F(2)=-Y(l)*Y(2)/YC3)+0. 140500  *CB/MOM  $  FC3)=Y(1) 
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FC4)=YC2)*YC5)/YC3)-YC6) 

FC5)=-2.0*YC2)*YC4)/YC3)+CYCl)*YC5))/YC3) 

FC6)=CCYC2)*YC2)-2.0*GM/YC3))*YC4)-YC1>K'YC2)*YC5)  )/ 

*CYC3)*YC3)) 

3  S=RKLDEQCN,Y,F,X,H,NT) 

IFCS-1. 0)79,6, 4 

IF  THE  FINAL  TIME  IS  UNKNOWN  IT  IS  USED  TO  STOP  THE  INTEGRATION. 

IF  THE  FINAL  TIME  IS  KNOW  IT  IS  US^D  HERE 
IF(X.QE.XI(3))5.2 
5  CONTINUE 

THE  WEIGHING  FACTION  IS  THE  SUM  OF  THE  SOUARES  OF  THE 
ERROR  IN  THE  FINAL  CONDITIONS.  HOWEVER  THE  SUM  OF  THE  ABSOLUTE  VALUE 
OF  THE  ERRORS  CAN  ALSO  BE  USED.  ALSO  THE  ERRORS  CAN  BE  WEIGHTED  IF 
THEY  DIFFER  IN  IMPORTANCE. 

A  WEIGHING  FUNCTION  OF  FINAL  CONDITIONS  MUST  BE  FORMULATED  HERE 
G=  YC1)*YC1)+CyC2)-0.8098)*CyC2)-0.8098)+CyC3)-1.525)* 

*(Y(3)-1 . 525) 

IFCY(1).LT. -20.0)  G=100 . 0 

1FCYCD.QT.20.0)  G-100.0  I 

79  CONTINUE 
PRINT  10. G 

10  FORMATCX1 7HOBUECTVNE  FUNC  =  ,  E20.10) 

DO  20  U=l, N 
20  PRINT  25,J,YCJ) 

25  FORMATCX3HYC  ,I2,3H)  =,  E20.10) 

DO  30  J=1.K 
30  PRINT  35  ,J,XICJ) 

35  FORMATCX4HXIC  ,I2,3H)  =,  E20.10) 

RETURN 
END 


FUNCTION  RKLDEOCN, Y,F, X,H,NT) 

RKLDQ 

1 

C 

D2  UCSD  RKLDEO 

F  63 

C 

MODIFIED  MAY  1963  CQ  REMOVED  FROM  CALLING  SEQUENCE) 

c 

TEST  OF  ALGOL  ALGORITHM 

DIMENSION  YC10)/FC10),OC10) 

RKLDQ 

2 

C 

REAL  X,H—  INTEGER  N,NT— COMMENT— BEGIN  INTEGER  1,J,L-REAL  A 
NT=NT+1 

RKLDO 

3 

GO  TO  Cl. 2, 3, 4), NT 

RKLDO 

4 

C 

GO  TO  SCNT) 

1 

DO  11  J=1,N 

RKLDQ 

5 

11 

QCU)=0. 

RKLDO 

6 

A= .  5 

RKLDO 

7 

X=X+H/2. 

RKLDO 

8 

GO  TO  5 

RKLDO 

9 
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2 

A=. 29289321881 

RKLDO  10 

GO  TO  5 

RKLDQ  11 

3 

A=l. 7071067812 

RKLDO  12 

X=X+H/2 . 

RKLDO  13 

GO  TO  5 

RKLDO  14 

4 

DO  41  1=1, N 

RKLDQ  15 

41 

YCl)=Y(0+H*FCO/6.-QCO/3. 

RKLDO  16 

NT=0 

RKLDO  17 

RKLDEO=2 . 

RKLDQ  18 

GO  TO  6 

RKLDO  19 

5 

DO  51  L=1,N 

RKLDO  20 

YCL)=YCL)+A*CH*FCL)-QCL)) 

RKLDQ  21 

51 

Q(L)=2.^*H^a)*<l.-3.*A)*QCL) 

RKLDO  22 

RKLDEQ=1. 

RKLDO  23 

6 

CONTINUE 

RKLDQ  24 

RETURN 

RKLDQ  25 

END 

RKLDQ  26 

[ 

*  PROGRAM  PENKT 

THIS  PROGRAM  WILL  SOLVE  A  TWO-POINT  BOUNDARY  VALUE  PROBLEM 
USING  A  DIRECT  SEARCH  SCHEME  (NELC  DIRECT)  AND 
THE  INTEGRATION  SCHEME  IS  THE  FOURTH  ORDER  RUNGE-KUTTA 
AS  MODIFIED  BY  GILL 

THE  XI  ARE  THE  INITIAL  CONDITIONS  THAT  ARE  UNKNOWN  AND  WHICH 
ARE  SPECIFIED  INITIALLY  IN  THE  PROGRAM  DRIVER 


K  REPRESENTS  THE  NUMBER  OF  UNKNOWN  VARIABLES. 

C  N  =  THE  NUMBER  OF  DIFFERENTIAL  EQUATIONS  TO  BE  INTEGRATED 
C 

C  THE  DELTA  REPRESENTS  A  STEP  CHANGE  IN  THE  INITIAL  CONDITIONS 
C  AND  WHICH  IS  USED  BY  THE  SEARCH  ROUTINE. 

C  THE  DELTA  SHOULD  BE  LARGE  INITIALLY.  THEN  IF  RESULTS  ARE  NOT 
C  SAT IS FACTORY .DELTA  CAN  BE  REDUCED. 

C 

C  THE  SUBROUTINE  FN  IS  USED  FOR  INTEGRATION  AND  ALSO  TO  COMPUTE 
C  A  WEIGHING  FUNCTION  OF  THE  FINAL  BOUNDARY  CONDITIONS 
C 

C  A  WEIGHING  FUNCTION  OF  FINAL  CONDITIONS  MUST  BE  FORMULATED 
C  THE  WEIGHING  FUNCTION  IS  USED  BY  THE  DIRECT  SEARCH  SCHEME 

C  TO  DETERMINE  THAT  THE  FINAL  CONDITIONS  ARE  MET 

C  THE  FINAL  CONDITIONS  ARE  WRITTEN  AS  EQUALITIES  P([Y)=0.0 

C  THE  WEIGHING  FUNCTION  IS  THE  SUM  OF  THE  SQUARES 

C  G=PCYCl))*P(Ya))+PCY(2))*P(Y(2))  +  ETC. 

C 

C  IF  THE  FINAL  TIME  IS  UNKNOWN,  AN  ADDITIONAL  VARIABLE  XI (K  )  is 

C  USED  TO  REPRESENT  THE  FINAL  TIME  AND 

C  INSERTED  IN  THE  INTEGRATION  STOPPING  CRITERION.  OTHERWISE 
C  THE  TRUE  FINAL  VALUE  OF  THE  INDEPENDENT  VARIABLE  IS  USED 
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c 

C  PROBLEMS  CAN  BE  ENCOUNTERED  WITH  STEP-SIZE  DELTAI  AND  THE 
C  INTEGRATION  STEP-SIZE  H.  THE  LOCAL  VS  GLOBAL  PROBLEM  IS  PRESENT 
C  OCCURRED  AND  IT  IS  NECESSARY  TO  RESTART  THE  ROUTINE 

C  WITH  NEW  INITIAL  CONDITIONS  TO  GUARANTEE  A  GLOBAL  OPTIMUM  SOLUTION. 
C 

C  ADVANTAGE  OF  THIS  ROUTINE  IS  THAT  IT  ONLY  REOUIRES  INTEGRATION 

C  OF  THE  DIFFERENTIAL  EQUATION  .  NO  PERTURBATION  EQUATIONS 

C  ARE  REQUIRED. 

C 

DIMENSION  EXC12),  XlCl2),OXa2),RXCl2),PX(12) 

EXTERNAL  RKLDEQ 

XI(1)=0.52 
XIC2)=0.3 
XI(3)=3 .03 

C  K=  NO.  OF  UNKNOWNS 

K=3 

C  INITIAL  DELTA 

DELTAS.  05 
DEL=0.001 
DELF=1 . 0E-06 
RHO=  0.625 

13  CALL  DIRECTCF  ,K,RHO,DEL,DELF,  XI,DELTA,EX,OX,RX,PX) 

END 
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r 
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SUBROUTINE  FN(G,K,XI) 

DIMENSION  XI  GO 
C 

C  THESE  QUANTITIES  ARE  NEEDED  FOR  INTEGRATION.  THE  ARRAYS  HAVE  THE 
C  SAME  ORDER  AS  THE  DIFFERENTIAL  EQUATIONS 
DIMENSION  Y(6),F(6) 

TYPE  REAL  MAGSQ,MAG/MOM 
C 

EXTERNAL  RKLDEQ 
N=6 
C 

C  THE  KNOWN  INITIAL  CONDITIONS  ARE  SPECIFIED  HERE 
YO)=0.0 
Y(2)=l .  0 
Y(3)=l . 0 
Y(4)=XIC1) 

Y(5)=XIC2) 

Y(6)=l. 0 
C 

C  THE  INITIAL  VALUE  OF  X,  THE  INDEPENDENT  VARIABLE  AND  THE  INTEGRATION 
C  STEP  SIZE  H  ARE  SPECIFIED  .  THE  RUNGE-KUTTA  INTEGRATION  SCHEME 
C  IS  SENSITIVE  TO  STEP  SIZE  AND  H  MUST  BE  CHOSEN  CAREFULLY 
H=0.030$  X=0.0 
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2  NT=0 

THE  RIGHT-HAND  SIDE  OF  THE  DIFFERENTIAL  EQUATIONS 
DY/DX  =  F(Y,X) 

IS  EVALUATED  NEXT 

6  mGSQ=YC4)*YC4)+Y(5;>*YC5) 

GM=1.0 

MAG=SQRTFCMAGSQ) 

MOM=l .0-0 . C749*X 
SB=  YC4)/MAG 
CB=  YC5)/MAG 

FCD=YC2)*YC2)/YC3)-GM/CYC3)#YC3))+0 . 140500  -*SB/MOM 
FC2)=-YC1)*YC29/YC39+0. 140500  *CB/MOM  $  F(3)=YCl) 
FC4)=YC2)*YC5)/YC3)-YC6) 

FC5)^-2.0*Y(2)*YC4)/YC3)+(YCl)5tlYC5))/YC3) 

FC6)=CCYC2)*Y(2)-2.0wGM/YC3))*YC4)-YCO*YC2)JfrYC5)  )/ 

*(Y(3)*YC30) 

C 

3  S=RKLDEQCN,Y,F,X,H,NT) 

IFCS-1. 0)79,6,4 
C 

C  IF  THE  FINAL  TIME  IS  UNKNOWN  IT  IS  USED  TO  STOP  THE  INTEGRATION. 

C  IF  THE  FINAL  TIME  IS  KNOW  IT  IS  USED  HERE 
4  1FCX.GE.XIC3))5,2 

5  CONTINUE 
C 

C  THE  WEIGHING  FUNCTION  IS  THE  SUM  OF  THE  SQUARES  OF  THE 
C  ERROR  IN  THE  FINAL  CONDITIONS.  HOWEVER  THE  SIM  OF  THE  ABSOLUTE  VALUE 
C  OF  THE  ERRORS  CAN  ALSO  BE  USED.  ALSO  THE  ERRORS  CAN  BE  WEIGHTED  IF 
C  THEY  DIFFER  IN  IMPORTANCE. 

C 

C  A  WEIGHING  FUNCTION  OF  FINAL  CONDITIONS  MUST  BE  FORMULATED  HERE 
G=  YCD*YCD+CYC2)-0.8098)*CyC2)-0.8098)+CYC3)-1.525)* 

*(YC3)-1. 525) 

IFCYCD.LT. -20. 09  G=100 . 0 
IFCYCD.QT.20.0)  G=100.0 
79  CONTINUE 
PRINT  10. G 

10  F0RNATCX17H08JECTVNE  FUNC  =  ,  E20.10) 

DO  20  J=1,N 

20  PRINT  25,J,YCU) 

25  FORMATCX3HYC  ,12,31-0  =,  E20.10) 

DO  30  J=1,K 
30  PRINT  35  ,J,XICJ) 

35  FORM4TCX4HXIC  ,I2,3H)  =,  E20.1Q) 

RETURN 

END 
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PROGRAM  TPBVPSCH 

THIS  PROGRAM  WILL  SOLVE  A  TWO -POINT  BOUNDARY  VALUE  PROBLEM 
USING  A  DIRECT  SEARCH  SCHEME  (UN I VAR)  AND 
THE  INTEGRATION  SCHEME  IS  THE  EXTRAPOLATION  METHOD  BY 
BOLESCH  AND  STOER 

THE  XI  ARE  THE  INITIAL  CONDITIONS  THAT  ARE  UNKNOWN  AND  WHICH 
ARE  SPECIFIED  INITIALLY  IN  THE  PROGRAM  DRIVER 

K  REPRESENTS  THE  NUMBER  OF  UNKNOWN  VARIABLES. 

N  =  THE  NUMBER  OF  DIFFERENTIAL  EQUATIONS  TO  BE  INTEGRATED 

THE  DELTAI  REPRESENTS  A  STEP  CHANGE  IN  THE  INITIAL  CONDITIONS 
AND  WHICH  IS  USED  BY  THE  SEARCH  ROUTINE. 

THE  DELTAI  SHOULD  BE  LARGE  INITIALLY.  THEN  IF  RESULTS  ARE  NOT 
SATISFACTORY, DELTAI  CAN  BE  REDUCED. 

THE  SUBROUTINE  FN  IS  USED  FOR  INTEGRATION  AND  ALSO  TO  COMPUTE 
A  WEIGHING  FUNCTION  OF  THE  FINAL  BOUNDARY  CONDITIONS 

A  WEIGHING  FUNCTION  OF  FINAL  CONDITIONS  MUST  BE  FORMULATED 
THE  WEIGHING  FUNCTION  IS  USED  BY  THE  DIRECT  SEARCH  SCHEME 
TO  DETERMINE  THAT  THE  FINAL  CONDITIONS  ARE  MET 
THE  FINAL  CONDITIONS  ARE  WRITTEN  AS  EQUALITIES  P(Y)=0.0 
THE  WEIGHING  FUNCTION  IS  THE  SUM  OF  THE  SQUARES 
G=P(Y(l))*PCYa))+PCY(2))*P(YC2))  +  ETC. 

IF  THE  FINAL  TIME  IS  UNKNOWN  ,  AN  ADDITION  VARIABLE  XI(K  )  IS 
USED  TO  REPRESENT  THE  FINAL  TIME  AND 

INSERTED  IN  THE  INTEGRATION  STOPPING  CRITERION.  OTHERWISE 
THE  TRUE  FINAL  VALUE  OF  THE  INDEPENDENT  VARIABLE  IS  USED 

PROBLEMS  CAN  BE  ENCOUNTERED  WITH  STEP-SIZE  DELTAI  AND  THE 
INTEGRATION  STEP-SIZE  H.  THE  LOCAL  VS  GLOBAL  PROBLEM  IS  PRESENT 
OCCURRED  AND  IT  IS  NECESSARY  TO  RESTART  THE  ROUTINE 

WITH  NEW  INITIAL  CONDITIONS  TO  GUARANTEE  A  GLOBAL  OPTiMUM  SOLUTION. 

ADVANTAGE  OF  THIS  ROUTINE  IS  THAT  IT  ONLY  REQUIRES  INTEGRATION 
OF  THE  DIFFERENTIAL  EQUATION  .  NO  PERTURBATION  EQUATIONS 
ARE  REQUIRED. 

DIMENSION  XI(4),X(4),DP(4),DEL(4)/XR(4) 

EXTERNAL  F 

XI C 1 5=  0.7 
XI (2)=0 . 3 
XI (3)=3 . 05 

K=  NO.  OF  UNKNOWS 
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INITIAL  DELTA 

DELTAI=0.04 

RH0= 1.618 

CALL  UN  I VAR(K, X I , RHO, DELTAI , X, DP , DEL, XR, FR) 
END 


SUBROUTINE  FNCG.K.XI) 

DIMENSION  XI OO 
EXTERNAL  F 

THESE  QUANTITIES  ARE  NEEDED  FOR  INTEGRATION,  THE  ARRAYS  HAVE  THE 
SAME  ORDER  AS  THE  DIFFERENTIAL  EQUATIONS 
DIMENSION  Y(6) 

TYPE  REAL  MAGSQ,MAG, MOM 

THE  KNOWN  INITIAL  CONDITIONS  ARE  SPECIFIED  HERE 
N=6 
YC19=0.0 
YC2)=1.0 
YC3)=1.0 
YC4)=XlCO 
YC5)=XIC2) 

YC6)=1.0 

THE  INITIAL  VALLE  OF  X,  THE  INDEPENDENT  VARIABLE  AND  THE  INTEGRATION 
STEP  SIZE  H  ARE  SPECIFIED  .  THE  RUNGE-KUTTA  INTEGRATION  SCHEME 
IS  SENSITIVE  TO  STEP  SIZE  AND  H  MUST  BE  CHOSEN  CAREFULLY 
H  IS  CHANGED  BY  THE  DIRECT  SEARCH  ROUTINE  HERE 
X=0.0$  H=XIC3)$ 


THE  RIGHT-HAND  SIDE  OF  THE  DIFFERENTIAL  EQUATIONS 
DY/DX  =  F(Y,X) 

IS  EVALUATED  NEXT 

1  CALL  DIFF  SYS(6,F,1.E-5,H,X,Y) 

C 

C  IF  THE  FINAL  TIME  IS  UNKNOWN  IT  IS  USED  TO  STOP  THE  INTEGRATION. 

C  IF  THE  FINAL  TIME  IS  KNOWN  IT  IS  USED  HERE 
5  CONTINUE 
C 

C  THE  WEIGHING  FUNCTION  IS  THE  SUM  OF  THE  SQUARES  OF  THE 
C  ERROR  IN  THE  FINAL  CONDITIONS.  HOWEVER  THE  SUM  OF  THE  ABSOLUTE  VALUE 
C  OF  THE  ERRORS  CAN  ALSO  BE  USED.  ALSO  THE  ERRORS  CAN  BE  WEIGHTED  IF 
C  THEY  DIFFER  IN  IMPORTANCE. 
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C 

C  A  WEIGHING  FUNCTION  OF  FINAL  CONDITIONS  MUST  BE  FORMJLATED  HERE 
G=  YCl)*YCl)+CYC2)-0.8098)«CYC2)-0.8098)+CY(3)-1.525)>!' 

*(Y(3)-1 . 525) 

79  CONTINUE 
PRINT  10, G 

10  FORMAT(X 1 7HOB JECTVNE  FUNC  =  ,  E20.10) 

DO  20  J=1,N 

20  PRINT  25, J, Y(J) 

25  FORMATCX3HYC  ,I2,3H)  =,  E20+10) 

DO  30  J-1,K 
30  PRINT  35  ,J,XlCJ) 

35  FORMAT(X4HXlC  ,I2,3H)  =,  E20.10) 

RETURN 


SUBROUTINE  F(X,Y,Z) 

TYPE  REAL  MAGSQ,  MAG, MOM 
DIMENSION  Y(6),Z(6) 

6  MAGS0=Y(9)frYC4)+YC5)i;cYC5) 

GM=  1 . 0 

MAG-SQRTFCMAGSQ) 

MOM=1.0-0.0749*X 
SE~  ;  C'O/MAG 
C8=  Y(5)/MAG 

Za)=Y(2)*Y(2)/YC3)^M/(Y(3)*YC3))+0. 140500  *SB/MOM 
ZC2)=-YCl)*YC2)/YC3)+0. 140500  *CB/MOM  $  ZC3)=YCl) 
Z(4)=Y(2)*Y(5)/YC3)-YC6) 

Z(5)=-2.0Y!YC2)*YC4)/YC3)+CYa)::‘Y(5))/YC3) 
ZC6)=CCYC2)^YC2)-2.0^GM,/Y(3))^YC4>YCD*YC2)"YC5)  )/ 

*(Y(3)*Y(3;» 

RETURN 

END 


REVERSE  SIDE  BLANK 


